set more off
set maxiter 100
capture log close
clear
cd "/Users/danielhill/Documents/ZiOP_jcr"
log using "ziopsims28Nov12-r.log", replace

/* program that generates data w/ dv following ziopr dgp (called in loop) */

capture program drop syndata
program syndata, rclass
	args obs b0 b1 b2 g0 g1 tau_1 tau_2
	clear
	set obs `obs'
	gen x1 = -5+(10)*runiform()
	gen x2 = rnormal(0,2)
	gen z1 = rbinomial(5,.5)
	scalar pee = .5
      matrix P = (1, pee \ pee, 1)
      drawnorm e1 e2, corr(P) 
	gen infl_lat = scalar(g0)+scalar(g1)*z1+e2
	gen infl = infl_lat > 0 
	sum infl, meanonly
	scalar _infl = r(mean)
	sum x1, meanonly
	scalar _meanx1 = r(mean)
	sum x2, meanonly
	scalar _meanx2 = r(mean)
	sum z1, meanonly
	scalar _meanz1 = r(mean)

	gen out_lat = scalar(b0)+scalar(b1)*x1+scalar(b2)*x2+e1
	gen outcome = out_lat > scalar(tau_1)
	replace outcome = 2 if out_lat > scalar(tau_2)

	gen y = infl*outcome
end

/* Create file and variables to store results from simulations */

postutil clear
postfile ziopsim obs sims b0 b1 b2 g0 g1 t2 r /// 
		     ob0 ob1 ob2 og1 ot2 ///
		     ovcv11 ovcv21 ovcv31 ovcv41 ovcv51 ///
		     ovcv22 ovcv32 ovcv42 ovcv52 ///
		     ovcv33 ovcv43 ovcv53 ///
		     ovcv44 ovcv54 ///
		     ovcv55 ///
		     zb0 zb1 zb2 zg0 zg1 zt2 ///
		     zvcv11 zvcv21 zvcv31 zvcv41 zvcv51 zvcv61 ///
		     zvcv22 zvcv32 zvcv42 zvcv52 zvcv62 ///
		     zvcv33 zvcv43 zvcv53 zvcv63 ///
		     zvcv44 zvcv54 zvcv64 ///
		     zvcv55 zvcv65 ///
		     zvcv66 ///
		     rb0 rb1 rb2 rg0 rg1 rt2 rho ///
		     rvcv11 rvcv21 rvcv31 rvcv41 rvcv51 rvcv61 rvcv71 ///
		     rvcv22 rvcv32 rvcv42 rvcv52 rvcv62 rvcv72 ///
		     rvcv33 rvcv43 rvcv53 rvcv63 rvcv73 ///
		     rvcv44 rvcv54 rvcv64 rvcv74 ///
		     rvcv55 rvcv65 rvcv75 ///
		     rvcv66 rvcv76 ///
		     rvcv77 ///
             meaninfl meanx1 meanx2 meanz1 rco rcz rcr ///
		     using "ziopsims28Nov12-r.dta", replace

/* Generate macros for number of observations in data, number of simulations, seed value*/

timer on 1
local obslist "4000"
local n = 5000
local sd = 666
local constant -6 -4.4 -3.3 -2.5 -1.7 -0.5 3 

/* Simulation loop */
	foreach obs of local obslist {
      local sd = `sd'
      qui set seed `sd'

      dis _n(2) as result "Number of Observations: `obs'"
      
	    foreach j of local constant {
      
			local i = 1

      		while `i' <= `n' {
	
			scalar b0 = 1
			scalar b1 = .5
			scalar b2 = 1.5
			scalar g1 = 1
			scalar g0 = `j'
			scalar tau_1 = 0
			scalar tau_2 = 2

		   	qui syndata `obs' scalar(b0) scalar(b1) scalar(b2) ///
        	scalar(g0) scalar(g1) /// 
			scalar(tau_1) scalar(tau_2)
	
			capture myoprob y x1 x2 z1
			matrix eb = e(b)
			matrix ev = e(V)
			local rco = _rc
			if `rco'==0 {
			scalar ob0 = eb[1,3] 
			scalar ob1 = eb[1,1]
			scalar ob2 = eb[1,2]
			scalar og1 = eb[1,4]
			scalar ot2 = eb[1,5]
			scalar ovcv11 = ev[1,1]
			scalar ovcv21 = ev[2,1]
			scalar ovcv31 = ev[3,1]
			scalar ovcv41 = ev[4,1]
			scalar ovcv51 = ev[5,1]
			scalar ovcv22 = ev[2,2]
			scalar ovcv32 = ev[3,2]
			scalar ovcv42 = ev[4,2]
			scalar ovcv52 = ev[5,2]
			scalar ovcv33 = ev[3,3]
			scalar ovcv43 = ev[4,3]
			scalar ovcv53 = ev[5,3]
			scalar ovcv44 = ev[4,4]
			scalar ovcv54 = ev[5,4]
			scalar ovcv55 = ev[5,5]
      		}
			else {
			scalar ob0 = . 
			scalar ob1 = .
			scalar ob2 = .
			scalar og1 = .
			scalar ot2 = .
			scalar ovcv11 = .
			scalar ovcv21 = .
			scalar ovcv31 = .
			scalar ovcv41 = .
			scalar ovcv51 = .
			scalar ovcv22 = .
			scalar ovcv32 = .
			scalar ovcv42 = .
			scalar ovcv52 = .
			scalar ovcv33 = .
			scalar ovcv43 = .
			scalar ovcv53 = .
			scalar ovcv44 = .
			scalar ovcv54 = .
			scalar ovcv55 = .
      		}
	
			capture ziop1 y x1 x2, infl(z1)
			matrix eb = e(b)
			matrix ev = e(V)
			local rcz = _rc
			if `rcz'==0 {
			scalar zb0 = eb[1,3] 
			scalar zb1 = eb[1,1]
			scalar zb2 = eb[1,2]
			scalar zg0 = eb[1,5]
			scalar zg1 = eb[1,4]
			scalar zt2 = eb[1,6]
			scalar zvcv11 = ev[1,1]
			scalar zvcv21 = ev[2,1]
			scalar zvcv31 = ev[3,1]
			scalar zvcv41 = ev[4,1]
			scalar zvcv51 = ev[5,1]
			scalar zvcv61 = ev[6,1]
			scalar zvcv22 = ev[2,2]
			scalar zvcv32 = ev[3,2]
			scalar zvcv42 = ev[4,2]
			scalar zvcv52 = ev[5,2]
			scalar zvcv62 = ev[6,2]
			scalar zvcv33 = ev[3,3]
			scalar zvcv43 = ev[4,3]
			scalar zvcv53 = ev[5,3]
			scalar zvcv63 = ev[6,3]
			scalar zvcv44 = ev[4,4]
			scalar zvcv54 = ev[5,4]
			scalar zvcv64 = ev[6,4]
			scalar zvcv55 = ev[5,5]
			scalar zvcv65 = ev[6,5]
			scalar zvcv66 = ev[6,6]
      		}
			else {
			scalar zb0 = . 
			scalar zb1 = .
			scalar zb2 = .
			scalar zg0 = .
			scalar zg1 = .
			scalar zt2 = .
			scalar zvcv11 = .
			scalar zvcv21 = .
			scalar zvcv31 = .
			scalar zvcv41 = .
			scalar zvcv51 = .
			scalar zvcv61 = .
			scalar zvcv22 = .
			scalar zvcv32 = .
			scalar zvcv42 = .
			scalar zvcv52 = .
			scalar zvcv62 = .
			scalar zvcv33 = .
			scalar zvcv43 = .
			scalar zvcv53 = .
			scalar zvcv63 = .
			scalar zvcv44 = .
			scalar zvcv54 = .
			scalar zvcv64 = .
			scalar zvcv55 = .
			scalar zvcv65 = .
			scalar zvcv66 = .
	        }

			capture ziop2 y x1 x2, infl(z1)
			matrix eb = e(b)
			matrix ev = e(V)
			local rcr = _rc
			if `rcr'==0 {
			scalar rb0 = eb[1,3] 
			scalar rb1 = eb[1,1]
			scalar rb2 = eb[1,2]
			scalar rg0 = eb[1,5]
			scalar rg1 = eb[1,4]
			scalar rt2 = eb[1,6]
			scalar rho = eb[1,7]
			scalar rvcv11 = ev[1,1]
			scalar rvcv21 = ev[2,1]
			scalar rvcv31 = ev[3,1]
			scalar rvcv41 = ev[4,1]
			scalar rvcv51 = ev[5,1]
			scalar rvcv61 = ev[6,1]
			scalar rvcv71 = ev[7,1]
			scalar rvcv22 = ev[2,2]
			scalar rvcv32 = ev[3,2]
			scalar rvcv42 = ev[4,2]
			scalar rvcv52 = ev[5,2]
			scalar rvcv62 = ev[6,2]
			scalar rvcv72 = ev[7,2]
			scalar rvcv33 = ev[3,3]
			scalar rvcv43 = ev[4,3]
			scalar rvcv53 = ev[5,3]
			scalar rvcv63 = ev[6,3]
			scalar rvcv73 = ev[7,3]
			scalar rvcv44 = ev[4,4]
			scalar rvcv54 = ev[5,4]
			scalar rvcv64 = ev[6,4]
			scalar rvcv74 = ev[7,4]
			scalar rvcv55 = ev[5,5]
			scalar rvcv65 = ev[6,5]
			scalar rvcv75 = ev[7,5]
			scalar rvcv66 = ev[6,6]
			scalar rvcv76 = ev[7,6]
			scalar rvcv77 = ev[7,7]
      		}
			else {
			scalar rb0 = . 
			scalar rb1 = .
			scalar rb2 = .
			scalar rg0 = .
			scalar rg1 = .
			scalar rt2 = .
			scalar rho = .
			scalar rvcv11 = .
			scalar rvcv21 = .
			scalar rvcv31 = .
			scalar rvcv41 = .
			scalar rvcv51 = .
			scalar rvcv61 = .
			scalar rvcv71 = .
			scalar rvcv22 = .
			scalar rvcv32 = .
			scalar rvcv42 = .
			scalar rvcv52 = .
			scalar rvcv62 = .
			scalar rvcv72 = .
			scalar rvcv33 = .
			scalar rvcv43 = .
			scalar rvcv53 = .
			scalar rvcv63 = .
			scalar rvcv73 = .
			scalar rvcv44 = .
			scalar rvcv54 = .
			scalar rvcv64 = .
			scalar rvcv74 = .
			scalar rvcv55 = .
			scalar rvcv65 = .
			scalar rvcv75 = .
			scalar rvcv66 = .
			scalar rvcv76 = .
			scalar rvcv77 = .
      		}

/* Post results */

			post ziopsim (`obs') (`n') (scalar(b0)) (scalar(b1)) (scalar(b2)) ///
			(scalar(g0)) (scalar(g1)) (scalar(tau_2)) (scalar(pee)) ///
            (scalar(ob0)) (scalar(ob1)) (scalar(ob2)) (scalar(og1)) (scalar(ot2)) ///
			(scalar(ovcv11)) (scalar(ovcv21)) (scalar(ovcv31)) (scalar(ovcv41)) ///
			(scalar(ovcv51)) ///
			(scalar(ovcv22)) (scalar(ovcv32)) (scalar(ovcv42)) (scalar(ovcv52)) ///
			(scalar(ovcv33)) (scalar(ovcv43)) (scalar(ovcv53)) ///
			(scalar(ovcv44)) (scalar(ovcv54)) ///
			(scalar(ovcv55)) ///
		    (scalar(zb0)) (scalar(zb1)) (scalar(zb2)) ///
			(scalar(zg0)) (scalar(zg1)) (scalar(zt2)) ///
			(scalar(zvcv11)) (scalar(zvcv21)) (scalar(zvcv31)) (scalar(zvcv41)) ///
			(scalar(zvcv51)) (scalar(zvcv61)) ///
			(scalar(zvcv22)) (scalar(zvcv32)) (scalar(zvcv42)) (scalar(zvcv52)) ///
			(scalar(zvcv62)) ///
			(scalar(zvcv33)) (scalar(zvcv43)) (scalar(zvcv53)) (scalar(zvcv63)) ///
			(scalar(zvcv44)) (scalar(zvcv54)) (scalar(zvcv64)) ///
			(scalar(zvcv55)) (scalar(zvcv65)) ///
			(scalar(zvcv66)) ///
		    (scalar(rb0)) (scalar(rb1)) (scalar(rb2)) /// 
	        (scalar(rg0)) (scalar(rg1)) (scalar(rt2)) (scalar(rho)) ///
			(scalar(rvcv11)) (scalar(rvcv21)) (scalar(rvcv31)) (scalar(rvcv41)) ///
			(scalar(rvcv51)) (scalar(rvcv61)) (scalar(rvcv71)) ///
			(scalar(rvcv22)) (scalar(rvcv32)) (scalar(rvcv42)) (scalar(rvcv52)) /// 
			(scalar(rvcv62)) (scalar(rvcv72)) ///
			(scalar(rvcv33)) (scalar(rvcv43)) (scalar(rvcv53)) (scalar(rvcv63)) ///
			(scalar(rvcv73)) ///
			(scalar(rvcv44)) (scalar(rvcv54)) (scalar(rvcv64)) (scalar(rvcv74)) ///
			(scalar(rvcv55)) (scalar(rvcv65)) (scalar(rvcv75)) ///
			(scalar(rvcv66)) (scalar(rvcv76)) ///
			(scalar(rvcv77)) ///
            (scalar(_infl)) (scalar(_meanx1)) (scalar(_meanx2)) (scalar(_meanz1)) ///
            (`rco') (`rcz') (`rcr')
	
			drop x1 x2 z1 e1 e2 infl_lat infl out_lat outcome y
			if mod(`i',1000)==0 display as result "`i' " _c
 
			local i = `i' + 1
  			}
		}
	}
postclose ziopsim	
set more on
log close

timer off 1
timer list 1
